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Abstract. Numerical simulations of the solar chromosphere have become increasingly 
realistic over the past 5 years. However, many observed chromospheric structures and be- 
havior are not reproduced. Current models do not show fibrils in Ca II 8542 A, and neither 
reproduce the Call 8542 A bisector. The emergent Ho- line core intensity computed from 
the models show granulation instead of chromospheric shocks or fibrils. I discuss these 
deficiencies and speculate about what physics should be included to alleviate these short- 
comings. 
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1. Introduction 

The solar chromosphere is in many ways the 
most complicated and most interesting region 
of the solar atmosphere. 

The chromosphere is the interface layer 
between the convection zone/photosphere and 
the corona where the value of the plasma 
y8-parameter changes from larger to smaller 
than unity. Consequently, its overall charac- 
ter changes from a convective overshoot layer 
into a regime where the three-dimensional 
(3D) magnetic topology determines the struc- 
ture. The chromosphere is convectively stable, 
so the dominant input of mechanical energy 
comes in the form of waves excited by gran- 
ules and the buffeting of magnetic concentra- 
tions in the photosphere. It also absorbs ra- 
diation from both the photosphere (e.g., H~, 
Balmer-continuum) and the corona (He con- 
tinua). Thermal conduction in the corona has 
a strong influence on the shape of the transi- 
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tion region, and hence of the upper chromo- 
sphere. The dominant chromospheric energy 
loss is through radiation in strong lines. Thus, 
a comprehensive model of the chromosphere 
cannot be made without modeling the under- 
lying photosphere and upper convection zone 
and the overlying lower corona. 

The wide range of pertinent physical pro- 
cesses require that a comprehensive descrip- 
tion of the chromosphere needs to incorporate 
a careful combination of atomic and molecular 
physics, radiative transfer, MHD and plasma 
physics. This complicated physics requires nu- 
merical modeling. 

From the above discussion I conclude that 
a "minimally satisfying self-consistent model" 
of the chromosphere should 

- be time-dependent; 

- be three-dimensional; 

- model radiative heating and cooling; 

- have a lower boundary in the convection 
zone, with a boundary condition that feeds 
enough energy into the simulation to main- 
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tain a the radiative energy loss at the solar 
effective temperature; 
- have an upper boundary in the corona. If 
the model cannot self-consistently main- 
tain a corona, it must provide enough en- 
ergy flux to keep plasma near the boundary 
at coronal temperatures. 

Within these constraints one can add features 
and physics to the model to increase its real- 
ism, while at the same time trying to keep the 
number of free parameters to a minimum. 

The above list does not imply that mod- 
els without those properties are somehow "bad 
models". To the contrary, simpler models are 
often instrumental in identifying and charac- 
terizing relevant physical processes. However, 
a thorough understanding of the chromosphere 
ultimately requires a self-consistent model that 
reproduces all observed chromospheric behav- 
ior. 



2. Current status 

The current state-of-the-art in self-consistent 
chromospheric modeling is the modeling 
performed with the radiation-magneto- 
hydrodynami cs (RMHD) Oslo Stagger 
Code (OSC iHansteenl 12004 and the MPI- 
parallelized Bifrost code curren tly under 
development (Gud iksen et al.ll20~ToT) . 

These codes are based o n the methods pio- 
neered bv lNordlundl (fl982) in 3D convection 
simulations. They use the single-fluid MHD 
approach to model the solar gas and mag- 
netic field. Radiative transfer in the photo- 
sphere and lower chromosphere is treated us- 
ing m ulti-group opa city binning with scatter- 
ing dSkartlienl |2000|) with typically 4 radi- 
ation bins. This multi-group scheme is not 
well suited to model the few strong chromo- 
spheric lines that cause the dominant radia- 
tive energy losses in the upper chromosphere. 
Therefore both codes employ parameterized 
radiative cooling in those lines in the up- 
per chromosphere based on the detailed radia- 
tive transfer in a ID dynamical chromospheric 
model (computed with the R ADYN code, see 
e.g., ICarlsson"& Stein 2002). Radiative losses 
in the corona are computed using an opti- 



cally thin radiative loss function. Heat conduc- 
tion parallel to the magnetic field is included. 
Two-dimensional simulations cannot maintain 
a corona by themselves and require addi- 
tional heating at the upper boundary. Three- 
dimensional models however self-consistently 
maintain a corona by Joule heating. 

Similar codes exist, such as Muram 
dVogler et al.ll2005h. C0 5 BOLD (iFrevtae et alJ 
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dAbbettll2008l) . However these codes so far do 
not support a corona (C0 5 BOLD, Muram) or 
lack detailed treatment of the radiation field 
(RADMHD). 



OSC has be en used to stu d y dyn amic 
fibrils (DFs) by iHansteen et al.l d2006l) and 
iDe Pontieu et alJ (l2007al) . They conclude that 
DFs are driven by magneto-acoustic shocks 
based on the similarity in dynamic behavior 
of wave-guided shocks in a 2D simulation and 
observations of DFs in the core of Ha. Their 
w ork has been extended to 3D in simulations 
bv lMartmez-Svkora et all (l2009bl) . 



3D simulations of flux emergence from 
the convection zone through the chromo- 
sphere into the corona were per formed by 
iMartmez-Svkora et alJ d2008l l2009al) A simi- 
lar simulation showed the presence of Alfvenic 
waves with properties very much like what 
has been observed with t he Hinode spacecraft 
dDe Pontieu et alj|2007bl) . 

For a good overview of 3D solar atmo- 
sphere modeling and post-processing in the 
form of 3D rad i ative transf er I refer the 
reader to lCarlsso n (2008, 2009) and references 
therein. 

Despite the impressive progress that has 
been made in the last 5 years, many, often sur- 
prisingly basic properties of the chromosphere 
are not yet reproduced in the models. In the 
next section I discuss some of those observa- 
tional facts and speculate about what is missing 
in the models. 
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3. Open questions and challenges 

3.1. Lack of extended magnetic 
structure in Ca II 8542 A 

The first thing one notices in Ca II 8542 A line- 
core images are the elongated structures de- 
lineating magnetic field lines around network 
regions. The NLTE Call 8542 A computation 
from a snapshot of a 3D s imulation of the chro- 
mosphere with OSC of iLeenaarts et al.l (I2009I) 
does, however, not show any such structures, 
despite the presence of a relatively strong mag- 
netic field. Instead it shows the signature of os- 
cillations and shocks, both in relatively field- 
free regions as in regions with stronger mag- 
netic fields. Similar structures are observed in 
Ca II 8542 A quiet regions away from the net- 
work. The lack of fibrils in the simulation is in 
my opinion the most pressing challenge mod- 
els of the chromosphere currently face. 

There can be many reasons for this lack. 
It might be that such fibrils require large-scale 
magnetic field configurations that do not fit in 
the typical size of the computational domain of 
16x16x16 Mm 3 . It might be that the simula- 
tion cannot form fibrils due to a lack of grid 
resolution. Alternatively, fibril-like structures 
might exist in the models, but the lack of reso- 
lution inhibits mass loading of the fibrils, pos- 
sibly due to too large numerical viscosity (see 
also Sec. 13.21 ). This lack of mass then trans- 
lates into optically thin fibrils. Another possi- 
bility is that Ca II 8542 A opacity in the fibrils 
is caused by physical processes currently not 
modeled in the simulations, for example non- 
equilibrium ionization or plasma effects that 
cannot be modeled within the MHD assump- 
tion. 

3.2. Line broadening and bisector 

The spatially averaged Ca II 8542 A line pro- 
file computed in 3D NLTE from the same OSC 
snapshot shows a too narrow and too deep 
line core compared to observati ons (see Figs. 2 
and 3 of Leenaarts et al.1 120091) . These calcu- 
lations did not employ any non-physical mi- 
croturbulence, and had a rather large grid cell 
size of 64 x 64 x 32 km 3 . The too narrow core 



was attributed to this lack of spatial resolu- 
tion. Indeed, a newer simulation by Carlsson 
(private communication) with a cell size of 
32 x 32 X 16 km 3 shows a shallower and wider 
average Ca II 8542 A core, but not yet as wide 
and shallow as observed. The increase in grid 
resolution causes an increase in amplitude of 
the velocity variations in the mid and upper 
chromosphere, and hence a widening of the av- 
erage profile. 

This result suggests that with even higher 
resolution the width of the line core might at 
least be qualitatively explained. However, the 
observed profile is at least partially formed in 
fibrils and as long as the models do not repro- 
duce these the core width and depth remains 
not completely understood. 

The low resolution simulation also did not 
reproduce the observed inverse C-shape of the 
Ca II 8542 A line bis ector. Simila r result s were 
obtained earlier by lUitenbroekl (l2006h who 
computed the NLTE line profiles from both 
a 3D simulation that did not extend into the 
corona and a dynamic ID model. It would be 
interesting to see whether the bisector shape 
can be reproduced by high-resolution simula- 
tions, or that it is determined by phenomena 
currently not present in the simulations. 

3.3. Lack of Ha opacity 

The identification of DFs as magneto-acoustic 
shocks guided along magnetic field lines has 
been established based on the excellent agree- 
ment of dynamical pro perties observed in Ha 
with 2D simulations (lHansteen et~ai1 120061: 
iDe Pontieu e"tai]|2007al) . 

In a simulation with identical reso- 
lution and magneti c field configuration, 
ILeenaarts et al.l d2007l) investigated the effect 
of non-equilibrium hydrogen ionization on 
the chromosphere. They found that the hy- 
drogen n = 2 level populations and column 
densities are much higher in the simulated 
DFs than elsewhere in the chromosphere. 
This is consistent with the observations where 
DFs appear dark on top of a brighter, deeper 
formed background. However, a subsequent 
NLTE column-by-column radiative transfer 
computation employing non-equilibrium 
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Fig. 1. Vertically emergent Ha line core intensity from a 2D radiation-MHD simulation with 
non-equilibrium hydrogen ionization taken into account. Startup effects are visible the first 3 
minutes. There is a weak imprint of upward propagating shocks, superimposed on a slower- 
evolving pattern caused by granulation and re versed granulation. No parabolic outlines of DFs 
are present as in the observations (see Fig. 7 of lDe Pontieu et ai1l2007al) . 
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Fig. 2. Results of a NLTE radiative transfer computation of the Ha line from a 3D radiation- 
MHD simulation with an LTE equation of state. The black line in panels a and b indicates the cut 
shown in panel c and d. Panel a: vertically emergent intensity in the far blue wing; b: vertically 
emergent intensity at line center; c: vertically emergent intensity along the cut indicated in panels 
a and b for the blue wing (grey, left-hand scale) and line center (black, right hand scale); d: t = 1 
height for the blue wing (grey) and line center (black). 



temperatures and electron densities, but as- 
suming statistical equilibrium for the hydrogen 
level-populations, showed that the upper chro- 
mosphere is transparent in the Ha line core in 
this simulation. 

Fig. Q] shows the time-evolution of the ver- 
tically emergent Ha line center intensity. There 
are magnetic field concentrations at x = 4 Mm 
and x = 12 Mm along which DFs run up in 



the chromosphere. Their tops should appear as 
parabolic curves in this xf-slice, but instead the 
slice shows a slowly evolving pattern of granu- 
lation and reverse granulation. Superimposed 
on this pattern is a weak imprint of upward 
propagating shocks (the diagonal curves ema- 
nating from the magnetic field concentrations, 
for example the one starting at (x, t) — (5.5, 10) 
Analysis of the contribution functions shows 
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that this imprint forms below 1 Mm, much 
lower than the typical height of DFs (see Fig. 
16 of lDe Pontieu et al.ll2007al) . 

Fig.|2]shows the results of a NLTE column- 
by-column radiative transfer computation of a 
sn apshot of a 3D simulation of fl ux-emergence 
by Ma rtfnez-Svkora et al.l d2008l) . This simula- 
tion used an LTE equation of state unlike the 
2D simulation that used a more realistic non- 
equilibrium equation of state. Panels a and b 
show that the Ha line core shows granulation, 
even though its r = 1 height lies between 0.3 
and 1.7 Mm (panel d). The line source func- 
tion is strongly scattering and set in the upper 
photosphere, where it is sensitive to the tem- 
perature pattern of the granulation. 

Both the 2D non-equilibrium simulation 
and the 3D LTE simulation show an Ha line 
center intensity pattern that is set in the upper 
photosphere, in stark contrast to the observa- 
tions which show fibrilar structure in all but the 
quietest areas of the atmosphere. Even in such 
quiet areas the observations do not show gran- 
ulation, but signatu res of Shockwaves instead 
dRutten et al.ll2008l) . 

Some of the reasons that might cause the 
absence of fibrils in the simulated Ca II 8542 A 
line core also apply here: too small compu- 
tational domain, too low spatial resolution or 
the limitations of the MHD assumption. The 
2D simulation necessarily had a rather simple 
magnetic field topology, which might inhibit 
the formation of fibrils and the mass loading 
of the upper chromosphere. The 3D simulation 
did not suffer from this, but it had lower reso- 
lution and did not model non-equilibrium ion- 
ization. High resolution 3D models with non- 
equilibrium hydrogen ionization might have a 
chromosphere that is optically thick in Ha. 

Hydrogen line formation is however noto- 
riously difficult to model. The non-equilibrium 
ionization of hydrogen does not only affect the 
electron density and temperature in the chro- 
mosphere, it also has a strong effect on hydro- 
gen line formation. The standard NLTE line 
formation assumption of statistical equilibrium 
is not valid because of the slow collisional and 
radiative transition rates as compared to the hy- 
drodynamical timescale. In addition at least the 
Lyman-a nd /3 lines need to be modeled with 



partial frequency redistribution (PRD) because 
of their strong influence on the Ha line. Both 
time-dependent radiative transfer and PRD ef- 
fects were ignored in the results of Figs. Q] 
and [2] Inclusion of these effects might signif- 
icantly increase the Ha opacity. 

Thus, proper 3D modeling of the Ha line 
requires full time-dependent radiative transfer 
with PRD in tandem with the MHD evolution, 
a Herculean task that has not even been done 
in ID hydrodynamic simulations. So far all 
ID simulations perform radiative transfer 
in co m plete redistribution ( Carlsson & Stein] 
| 2002t [Rammacher & Ulmschneiderl 2003t 
iKasparova et al.ll2009l) . 

An easier way out is to perform the 3D 
radiation-MHD simulations with approximate 
non-equilibrium hydrogen ionization follow- 
ing the method develo ped by ISolluml i 19991) 
as was done in 2D by iLeenaarts et al.1 (2007) 
in order to get reasonably accurate values for 
the electron density and temperature. One can 
then perform time-dependent radiative radia- 
tive transfer with PRD on a series of snapshots 
from the simulation to get the Ha line pro- 
files. This separates the detailed radiative trans- 
fer from the MHD computation, which signif- 
icantly reduces the difficulty of the problem. 
However, even with this simplification it re- 
mains a daunting task to accurately model Ha 
line formation in the dynamic chromosphere. 

4. Summary & conclusions 

Multidimensional radiation-MHD simulations 
of the solar chromosphere are becoming in- 
creasingly realistic. They have, amongst other 
things, been used to study dynamic fibrils, flux 
emergence and the propagation of Alfvenic 
waves into the corona. 

The models yet fail to reproduce some ob- 
servational facts: 

- simulated Ca II 8542 A core images do not 
show fibrils; 

- the average simulated Ca II 8542 A core is 
deeper and narrower than observed and its 
bisector not show the observed inverse C- 
shape; 

- the simulated Ha line core shows granu- 
lation, both in 2D simulations with non- 
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equilibrium effects taken into accounts and 
in 3D with an LTE equation of state. 

In order to reproduce these the models 
probably require higher resolution, larger com- 
putational domains, an improved treatment of 
radiation and non-equilibrium hydrogen ion- 
ization and inclusion of plasma effects. In or- 
der to model hydrogen transitions, in particular 
Ha, 3D NLTE time-dependent radiative trans- 
fer codes including PRD need to be developed. 

Acknowledgements. J.L. acknowledges financial 
support by the European Commission through 
the SOLAIRE Network (MTRN-CT-2006-035484) 
and the Netherlands Organisation for Scientific 
Research (NWO). 



References 

Abbett, W. P. 2008, in Astronomical Society 
of the Pacific Conference Series, Vol. 383, 
Subsurface and Atmospheric Influences on 
Solar Activity, ed. R. Howe, R. W. Komm, 
K. S. Balasubramaniam, & G. J. D. Petrie , 
327 

Carlsson, M. 2008, Physica Scripta Volume T, 

133,014012 
Carlsson, M. 2009, Memorie della Societa 

Astronomica Italiana, 80, 606 
Carlsson, M. & Stein, R. F. 2002, ApJ, 572, 

626 

De Pontieu, B., Hansteen, V. H., Rouppe van 

der Voort, L., van Noort, M., & Carlsson, M. 

2007a, ApJ, 655, 624 
De Pontieu, B., Mcintosh, S. W., Carlsson, M., 

et al. 2007b, Science, 318, 1574 
Freytag, B., Steffen, M., & Dorch, B. 2002, 

Astronomische Nachrichten, 323, 213 
Gudiksen, B. V., Carlsson, M., Hansteen, V., 

et al. 2010, A&A, in preparation 
Hansteen, V. H. 2004, in IAU Symposium, ed. 

A. V. Stepanov, E. E. Benevolenskaya, & 

A. G. Kosovichev, 385-386 
Hansteen, V. H., De Pontieu, B., Rouppe van 

der Voort, L., van Noort, M., & Carlsson, M. 

2006, ApJ, 647, L73 
Kasparova, J., Varady, M., Heinzel, P., 

Karlicky, M., & Moravec, Z. 2009, A&A, 

499, 923 



Leenaarts, J., Carlsson, M., Hansteen, V, & 
Rouppe van der Voort, L. 2009, ApJ, 694, 
L128 

Leenaarts, J., Carlsson, M., Hansteen, V, & 

Rutten, R. J. 2007, A&A, 473, 625 
Martmez-Sykora, J., Hansteen, V, & Carlsson, 

M. 2008, ApJ, 679, 871 
Martmez-Sykora, J., Hansteen, V, & Carlsson, 

M. 2009a, ApJ, 702, 129 
Martmez-Sykora, J., Hansteen, V, DePontieu, 

B., & Carlsson, M. 2009b, ApJ, 701, 1569 
Nordlund, A. 1982, A&A, 107, 1 
Rammacher, W. & Ulmschneider, P. 2003, 

ApJ, 589, 988 
Rutten, R. J., van Veelen, B., & Siitterlin, P. 

2008, Sol. Phys., 28 
Skartlien, R. 2000, ApJ, 536, 465 
Solium, E. 1999, Master's thesis, University of 

Oslo 

Uitenbroek, H. 2006, ApJ, 639, 516 

Vogler, A., Shelyag, S., Schiissler, M., et al. 

2005, A&A, 429, 335 
Wedemeyer, S., Freytag, B., Steffen, M., 

Ludwig, H.-G., & Holweger, H. 2004, 

A&A, 414, 1121 



